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We explore a signature of phase correlations in Fourier modes of dark mat- 
ter density fields induced by nonlinear gravitational clustering. We compute the 



distribution function of the phase sum of the Fourier modes, 9 kl + 6 k2 + ks , for 

■ triangle wavevectors satisfying fci + &2 + £3 = 0, and compare with the analytic 

prediction in perturbation theory recently derived by one of us. Using a series 

00 ■ of cosmological iV-body simulations, we extensively examine the time evolution 

O 

m 

Overall we find that the numerical results are remarkably consistent with the 
analytic formula from the perturbation theory. Interestingly the validity of the 
perturbation theory at a scale k corresponding to the wavevector k is determined 
by P{k)/V saimp , the ratio of the power spectrum P(k) and the sampling volume 



V^amp, not by k 3 P(k) as in the case of the conventional cosmological perturbation 
theory. Consequently this statistics of phase correlations is sensitive to the size of 



the sampling volume itself. This feature does not show up in more conventional 
cosmological statistics including the one-point density distribution function and 
the two-point correlation functions except as a sample-to-sample variation. Sim- 
ilarly if the sampling volume size V sajmp is fixed, the stronger phase correlation 
emerges first at the wavevector where P{k) becomes largest, i.e., in linear regimes 
according to the standard cosmological perturbation theory, while the distribution 
of the phase sum stays fairly uniform in nonlinear regimes. The above feature 
can be naturally understood from the corresponding density structures in real 
space as we discuss in detail. 
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1. Introduction 

The most conventional statistics in cosmology include the two-point correlation func- 
tion in real space and the power spectrum in Fourier space. They have been investigated 
both by numerous analytical and numerical methods and provided useful insights into the 
cosmological parameters and properties of the galaxy biasing among others mainly from the 
large-scale structure of the universe. The present universe, however, has much more com- 
plex structures, such as filaments, voids, and superclusters, which cannot be fully described 
by the two-point statistics. Therefore other statistics beyond the two-point statistics are 
desperately required for better understanding of the universe. 

The Fourier transform fj, of the density fluctuations, S(x) = (p(x) — p)/p, where p{x) 
is the density field and p is the spatial mean of the density field, is expressed in terms of the 
modulus |/fc| and the phase 6^ as follows: 

fk = |/fe|exp(i0 fc ). (1) 

The two-point statistics is defined in terms of the modulus alone, and thus the statistics of 
phases play complementary roles in characterizing the nature of the density field. Especially, 
the correlation of phases among different Fourier modes is a key ingredient in understanding 
the emergence of the non-Gaussian signature from the primordial Gaussian density field. 

In fact several statistics which carry the phase information have been already pro- 
posed in cosmology, including the void probability function (White 1979), the genus statistics 
(Gott et al. 1986), the Minkowski functionals (Mecke et al. 1994; Schmalzing & Buchert 1997). 
However, finding useful statistics of the Fourier phase itself is difficult mainly due to the cyclic 
property of the phase. For example, the one-point phase distribution turns out to be es- 
sentially uniform even in a strongly non-Gaussian field (Suginohara & Suto 1991) and then 
one cannot extract any useful information out of it. For this reason, previous studies of 
the Fourier phase have been mainly devoted to the evolution of phase shifts in individual 
modes (Ryden & Gramann 1991; Soda & Suto 1992; Jain & Bertschinger 1998)), and the 
phase differences between the Fourier modes (Scherrer et al. 1991; Coles & Chiang 2000; 
Chiang 2001; Chiang et al. 2002; Watts et al. 2003). There is still, however, a very incom- 
plete understanding of how phase correlations among different modes start to show up, or 
for the corresponding structure in real space that the strong phase correlation indicates. 

The connection between the higher-order statistics and the phase correlations has been 
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also suggested (Bertschinger 1992; Watts & Coles 2003). Recently, one of the present au- 
thors obtained an analytic expression for the distribution function of the "phase sum" 
#fci + #fc 2 + ' ' ' + ®k N where the corresponding wavevectors satisfy ki + k 2 + ■ ■ ■ + k N = 
(Matsubara 2003b). He discovered the general relation between the distribution of the phase 
sum and the hierarchy of polyspectra in the perturbation theory. Following his analytic 
results, we extensively study the behavior of the distribution of the phase sum of triangle 
wavevectors using a series of cosmological iV-body simulation. 

The plan of this paper is as follows. In §2, we briefly review the analytic expansion 
formula of the phase correlations with particular emphasis to the sampling volume depen- 
dence on the distribution function of the phase sum. Section 3 summarizes the simulation 
data and the computation method, and the results are shown in §4. Finally §5 is devoted to 
conclusions and discussion. 



2. Series expansion of the Distribution of the Phase Sum 

Starting from the Edgeworth-like expansion of the joint probability distribution function 
(PDF) (Matsubara 1995; Matsubara 2003a) of arbitrary sets of the modulus factor \f k \ and 
the phase factor 6 k in the Fourier mode fc, Matsubara (2003b) derived the joint PDF of phases 
among different Fourier modes of arbitrary closed wavevectors after integrating over the other 
variables. The condition of the closed wavevectors ensures the translational invariance of the 
statistics. In the lowest order approximation, the joint PDF of the three phases in triangle 
wavevectors kl ,9 k2 an d #fc 3 = —9 kl+k2 (&3 — ~~ &i ~~ &2, &i 7^ ^2) i n a sampling volume 
Kamp is written (Matsubara 2003b) as 

71-3/2 

|V S amp) OC 1 + ^j— p (3) (fci, k 2 \V samp ) cos(^fc 1 + k2 — 9 kl+k2 ), (2) 

P i3) (k 1} k 2 \V samp ) = B(k u k 2 ) 

P VV aamp P(|fe 1 |)P(|fe 2 |)P(|fc 1 + fc 2 |) 
For any periodic distributions on [0, 2ir], it is natural to expect that its lowest-order expansion 
gives a constant plus the cosine term. Therefore the most important point in the above 
formula lies in the fact that the proportional factor is given by the cumulant defined in 
terms of the bispectrum B(ki,k 2 ) and the power spectrum P(k) as follows: 

(/*/*> = % + k'P(\k\), (4) 
(/ fel / fe2 / fe3 ) = V^gS^+^Bikukt), (5) 

where (• • •) is the ensemble average and 8 k denotes the Kronecker delta which is 1 only 
when k = 0. We use the convention of the Fourier transform f k = V~ l l 2 J v d 3 xe~ lk ' x f(x) 
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in a finite box-size V. We explicitly write the volume dependence of p^\ which is derived 
in Appendix A. Equation (2) implies that the PDF of the phase sum, 9 kl + 9 k , 2 — 9 kl+k2 , is 
written as 

71-3/2 

V(6 kl + 6 k2 - 6> fel+fe2 |y samp ) oc 1 + -^-p (3) (fci, fc 2 |\/ samp ) cos(6» fcl + 6 k2 - 9 kl+k2 ). (6) 

We assume the statistical isotropy in the field distribution in the following analysis and 
thereby the distribution of the phase sum depends on \k±\, \k 2 \ and 6*12 in a fixed box-size 

^samp- 

If the hierarchical clustering ansatz is valid, B(ki,k 2 ) ~ P(ki)P{k 2 ) and thereby p^ 
is approximately given by >JP(k)/V samp . Therefore, the lowest-order approximation (eq. [6]) 
breaks down at a scale k satisfying the condition P(/c)/Vs amp > 1, which is different from 
the condition of the nonlinear clustering, k 3 P(k) > 1. When V^mp is sufficiently large, the 
higher-order terms p( N ^ become negligible and equation (6) is applicable even in nonlinear 
regimes in the sense of the conventional cosmological perturbation theory. 

The unusual feature of the statistics of the phase sum is that the degree of the non- 
uniformity explicitly depends on the size of the sampling volume Vs amp . This comes from 
the fact that the Ath-order terms p( N ^ have different volume dependence, oc Vsamjf 7 ' 2 (see 
Appendix A). Therefore, the distribution of phases in an infinite volume is always random 
in a statistical sense. The concept of the phase distribution function is meaningful only when 
averaged over many subsamples of a fixed Vs amp within the whole sample. 



3. A-body Simulation Data and Computational Method 

In order to study the behavior of the distribution of the phase sum, we apply the analytic 
formula (eq. [6]) on dark matter distribution using a series of P 3 M A-body simulations (Jing 
1998; Jing & Suto 1998). These simulations employ 256 3 particles in a cubic box with 256 3 
meshes and start from Gaussian-random distribution with periodic boundary conditions. We 
adopt a variety of the initial power spectra including four scale-free models with power-law 
indices n = —2, —1, 0, and +1. The simulations were terminated at the scale factor a = 1.0 
when the r.m.s. mass fluctuation smoothed over the top-hat filter with the scale of one-tenth 
of a box-length reached unity. We use the scale-free models at different scale factors so as 
to study the evolution of the distribution of the phase sum. We also analyze three cold 
dark matter (CDM) models, Lambda CDM (LCDM), Standard CDM (SCDM), and Open 
CDM(OCDM) at z — 0. The model parameters of all above simulations are summarized in 
Table 1. 
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Table 1. Simulation model parameters: the dimensionless matter-density parameter tt ; 
the dimensionless cosmological constant A ; the shape parameter T of the CDM transfer 
function (Bardeen et al. 1986); the r.m.s. density fluctuation amplitude smoothed by 
top-hat filter with the scale of 8/i _1 Mpc, a 8 . 



Model 


^0 


Ao 


r 


0"8 


LCDM 


0.3 


0.7 


0.21 


1 


SCDM 


1 





0.5 


0.6 


OCDM 


0.3 





0.25 


1 


scale-free (n = —2, —1, 0, +1) 


1 










Table 2. The original simulation box-size Vb ox and the ratio of sampling volume to 

original box-size, V samp /V hox . 



Model 




^samp / ^box 


CDM (LCDM,SCDM,OCDM) 
CDM (LCDM,SCDM,OCDM) 
scale-free (n = —2, —1, 0, +1) 


(SOO/i^Mpc) 3 
(lOO/i^Mpc) 3 


1 

1, 1/8, 1/64 
1, 1/8, 1/64 
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We analyze three different realizations for almost all of scale-free models and CDM 
models to estimate the sample-to-sample variance. In order to study the box-size dependence 
of phase correlations, we consider the CDM data with different box-size Vb ox of (300/i -1 Mpc) 3 
and (lOO/i^Mpc) 3 . Furthermore, we divide the original box data based on all of scale-free 
models and CDM models with volume of (lOO/i^Mpc) 3 into several sub-boxes with an equal 
volume V saxnp : 8 cubic boxes with a half length, and 64 cubic boxes with a quarter length of 
the original box, as shown in Table 2. 

The computation of the PDF of the phase sum proceeds as follows: first we construct 
the density fields defined on grid points using the cloud-in-cell interpolation from particle 
distribution, and then Fourier-transform them into /c-space. Next we choose two wavevectors 
ki and k 2 whose absolute values and open angle are within a certain range. In what follows, 
we consider, for simplicity and definiteness, only the configuration of |fei| and |fc 2 | whose 
absolute values are within the same range. Finally we obtain three phases of the Fourier 
modes for all sets of the chosen wavevectors, ki, k 2 and ki + k 2 and then compute the 
distribution of the phase sum 9^ + 0f~ 2 — Ok 1 +k 2 divided by the mean number per binned 
phase sum. We set a bin width of O.ln in 2ir range. We estimate in two ways; one is 
based on a fit of the PDF of the phase sum to the analytic formula of equation (6), and the 
other uses the direct definition (eq. [3]) in terms of the power spectrum and the bispectrum 
averaged over the same range of |fci|, \k 2 \ and 9i 2 . 

Throughout the paper the modulus of each wavevector is expressed in units of the 
Nyquist wavenumber of the grid in each field unless the units of the scale are explicitly men- 
tioned. For instance, k = [a, b] implies that the two wavevectors satisfy a < |fci|//cN y quist < b, 
and a < |fc 2 |/A:Nyquist < b- The Nyquist wavenumber is defined by 

, - 27T Tl g r id . . 

^Nyquist — — JJ^ Z~ , {< ) 

'samp 

where n gri d is the number of grids in a length of each box and Vj~Jm P is the length of the 
sampling box. In scale-free models, we choose n gri d = 256 x (V^amp/Vbox) 1 ^ 3 s° as to fix 
^Nyquist as 128 x 2tt/V^, where Vt, ox is the original simulation volume. We also confirm that 
the choice of n gr id does not change the results which we present in the next section. 
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4. Results 

4.1. Dependence of the phase sum PDF on the nonlinear gravitational 
evolution, scale and angle of the wavevectors, and the sampling volume 

size 

Figure 1 shows the time evolution, and the dependence on k, Q\2 and V saxav , of the 
distribution function of the phase sum (symbols) for scale-free models with power-law indices 
n = —2, —1, 0, and +1. For comparison, we plot with lines the lowest-order analytic formula 
(eq.[6]) with p^ directly evaluated from power spectra and bi-spectra of simulations. We 
note that the Poisson error of the distribution of the number of binned phase sum is very 
small. 

The symbols in Figure 2 plot the cumulant p^ as a function of the scale factor a, 
the central value of the wavevector k c (k = [k c — 0.05, k c + 0.05]), 9±2, and V^ amp . The 
corresponding lines indicate p^ evaluated from a fit of the phase-sum distribution to the 
analytic formula of equation (6). We choose the fitting range as 6^ +0k 2 —0k 1 +k 2 — [1-97T, 2n] 
and [On, O.ln] where the deviation from the lowest-order approximation is clearly seen. The 
quoted error-bars indicate the sample variance for n = — 2 model. 

Figures 1 and 2 employ a fiducial set of parameters, a — 1, k — [0.4,0.5], 9u = 
[120°, 130°], and = Vb ox , except for one parameter that is examined in each panel. 

The top panels in Figure 1 show the evolution of the PDF of the phase sum. In the 
earliest epochs, the density fields are still close to Gaussian and the PDF of the phase 
sum remains almost uniform. While the non-uniformity of the PDF develops at subsequent 
epochs, the upper-left panel in Figure 2 indicates that the lowest-order perturbation result 
well agrees with the simulation data even up to a = 1. 

The second-upper panels in Figure 1 show the scale dependence of the phase-sum PDF. 
Again the overall agreement between the lowest-order perturbation theory and simulation 
data is very good. Nevertheless only on large scales (k = [0.0, 0.1]), the lowest-order approx- 
imation starts to become invalid (see also the upper-right panel of Fig. 2). The fluctuations 
on large scales should be closer to the linear stage, and thus this behavior seems to be in- 
consistent with a naive expectation that non-Gaussianity should show up first from small 
scales. This apparent discrepancy can be reconciled by the fact that the cumulant p^ is 
approximately proportional to y/ P(k) /V samp as we emphasized in §2. In all the current data 
of the scale-free models, the nonlinearly evolved power spectra P(k) decrease with k (Fig. 
3). Therefore, as long as V^p is fixed, p^ is smaller for larger k, and thus the lowest-order 
perturbation provides a better approximation. This is an interesting and unique feature of 
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Fig. 1. — Various parameter dependence of the distribution of the phase sum 9k 1 +8k 2 ~ #fci+fc 2 
plotted by symbols in comparison with the perturbative formula (eq.[6]) plotted by lines for 
scale- free models. The power-law index n decreases from left to right. Each panel shows the 
time evolution, the scale dependence, the 6*12 dependence, and the sampling volume V samp 
dependence from top to bottom. 
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Fig. 2. — Comparison of p^ 3 \ki, ^2) plotted by symbols with the amplitude of the distribu- 
tion of phase sum at 9^ + 9k 2 — &ki+k 2 — divided by the constant factor 7r 3 / 2 /4 (see eq.[6]) 
plotted by lines for scale free models. The error-bars indicate the sample variances of 
and the amplitude of the distribution of phase sum for n = —2 model. Each panel shows the 
dependence of scale factor a (Upper-Left), scale( Upper- Right) , 9 12 (Lower-Left) and sampling 
volume Vsamp (Lower-Right). 
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Fig. 3. — Power spectrum P(k) in units of Vb ox & t various epochs in scale free models with 
n = 1( Upper- left), n = O(Upper-right), n = —l(Lower-left) and n = —2(Lower-right) . 
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the phase sum statistics that we examine here. 

The second-lower panels in Figure 1 show the 9 12 dependence of the phase-sum distribu- 
tion. As indicated more clearly in the lower-left panel of Figure 2, the value of p( 3 \ and thus 
the non-uniformity of the phase-sum PDF, increases slightly at 9± 2 > 120° for n = — 1 and 
—2 models. We suspect that this is ascribed to the presence of the large-scale filamentary 
structures in those models which the large 9i 2 (therefore small \ki + k 2 \) in /c-space may 
correspond to. 

The bottom panels in Figure 1 show the dependence on the sampling volume size V samp . 
Again unlike the conventional statistics, the PDF of the phase sum is indeed sensitive to the 
choice of K; am p. The lower-right panel of Figure 2 illustrates that p^ ^V ss , mp is approximately 
constant implying the validity of the hierarchical clustering ansatz as discussed in §2. This is 
why non-Gaussian feature of the PDF becomes substantially stronger and the higher-order 
terms become more important as V^ amp becomes smaller. 

Figures 4 and 5 plot the same results as Figures 1 and 2 but for the CDM models. The 
distribution of the phase sum is more uniform in SCDM models compared with other models 
mainly due to the smaller value of og (Table 1). This is because the small <jg corresponds to 
the small power P{k) with all over the scales and then p^ ~ yJP(k)/V has low value. The 
behavior of the parameter dependence is almost same as that in scale-free models. When a 
length of the simulation box is 300/i _1 Mpc, the lowest order approximation well agrees with 
simulations for any CDM models, while the approximation becomes worse as the box length 
decreases. 



Fig. 4. — Comparison between the distribution of phase sum 6 kl + 6 k2 — kl+k2 plotted 
by symbols and the lowest order perturbative equation (eq.[6]) plotted by lines for CDM 
models. The sampling volume V^mp is (300/i _1 Mpc) 3 , (lOO/i^Mpc) 3 , (50/i _1 Mpc) 3 , and 
(25/i _1 Mpc) 3 from top to bottom. The range of the modulus \k±\ and \k 2 \ is [0.1, 0.2] (open 
circles and solid lines), [0.4, 0.5] (closed triangles and dotted lines), and [0.8, 0.9] (crosses and 
dashed lines) in units of the Nyquist wavenumber /cNyquist written in each panel. 



4.2. Correspondence between phase correlations and density peaks 

As we emphasized over and over again, the PDF of the phase sum does depend on the 
size of the sampling volume. In order to understand the physical correspondence between 
the real-space distribution and the phase correlation, we compare the particle distribution 
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the dependence of V^ amp (Left) and |fei|(|fe 2 |) m units of Nyquist = 27r/(0.78/?,~ 1 Mpc) (Right). 
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from small sampling volumes on a sample-to-sample basis. 

First note that the non- uniformity of the phase sum arises from a few dominating density 
peaks in the sampling volume. Consider a field which has just one spherical peak with a 
density profile h(r) located at x . The density field in this idealized system is given by 
p(x) = h(\x — x \), and irrespective of the density profile of the peak, all the Fourier phases 
are aligned to 9k = — k ■ x . Thus the phase sum becomes zero in any combination of closed 
wavevectors. If another peak exists at a different position, the alignment of the phases is 
broken depending on the relative strength of the two peaks. If the peak height of the second 
peak is much smaller than that of the first, the phase alignment is almost preserved and the 
phase sum remains distributed around zero. On the other hand, if they have a comparable 
peak height, the phase alignment is significantly broken. In reality, all density peaks should 
contain substructure to some extent and many other structures in the field also contribute 
the eventual phase distribution. Nevertheless the above qualitative picture proves to be 
useful in understanding the physical meaning of the phase sum distribution. 

Motivated by the above consideration, we examine in detail the connection between 
the particle distribution and the PDF of the phase sum. Figure 6 plots the projected par- 
ticle distribution and the corresponding PDF of the phase sum for 8 different samples of 
(25/i _1 Mpc) 3 box from one of the LCDM simulations of (100/i _1 Mpc) 3 box. 

After smoothing the particle distribution over the Gaussian window of l/i _1 Mpc radius, 
we compute genus using CONTOUR 3D (Weinberg 1988). Then we define the peak heights 
of the 1st, 2nd and 3rd peaks, Si, 5 2 , and 5 3 by the threshold density contrasts which the 
number of genus in the box changes from from —1 to 0, from —2 to —1 and from —3 to 
—2, respectively. Those values are indicated in each panel of Figure 6. As expected from 
the above argument, Figure 6 clearly shows that the degree of non-uniformity of the PDF 
is related to the presence of the dominating peak in the box; the sample (d) has a single 
prominent density peak, and the phase sum is strongly distributed around zero, while the 
sample (b) has three almost comparable peaks resulting a fairly uniform PDF of the phase 
sum. 

To be more quantitative, we compute p^ and the ratio of the peak heights of the first 
and the second peaks for 64 samples of V^mp = (25/i _1 Mpc) 3 and 8 samples of V sa , mp = 
(50/i _1 Mpc) 3 , both constructed from one realization of LCDM model. The result plotted 
in Figure 7 indeed supports the strong correlation between the value of and the ratio 
of the peak heights. In fact, this naturally explains the dependence of the PDF on V samp . 
As K, amp increases, the number of density peaks in the sampling volume increases and those 
peaks tend to weaken the correlation of the phase sum. This clearly illustrated in Figure 8 
which plots the result for the box of K, amp = (50/i _1 Mpc) 3 comprising all the eight samples 
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of V^amp = (25/i _1 Mpc) 3 plotted in Figure 6. 



Fig. 6. — Comparison of the distribution of the phase sum with the value of three highest 
density peaks for eight fields with a box-length of 25/i~ 1 Mpc. The projected maps of dark 
matter particles are shown together. The distributions of the phase sum (symbols) are 
plotted together with the corresponding lowest order perturbative formula (lines) at three 
ranges of scale |fei| and |&2| i n [0.1,0.2] (filled circles, solid lines), [0.4, 0.5] (open triangles, 
dotted lines) and [0.7,0.8] (crosses, dashed lines) in units of /c Nyquist = 27r/(0.78/i~ 1 Mpc). 



5. Conclusions and Discussion 

We have performed a detailed numerical analysis of the phase correlation induced by 
nonlinear gravity in the universe. Following recent analytic work by Matsubara (2003b) 
in perturbation theory, we have computed the distribution function of the phase sum in 
triangle wavevectors and compared with his analytic formula. Using a large set of iV-body 
simulations, we have explored the behavior of phase correlations for various configurations of 
wavevector triangles in the Fourier space at different epochs in a range of different sampling 
volumes. We find that the agreement of the simulations and the analytic perturbation 
formula is generally good. 

Quite interestingly the series expansion of the analytic formula breaks down for P(k) /V sg , n 
1, which is in marked contrast to the conventional linear /nonlinear criterion, P(k)k 3 ~ 1. 
Therefore the statistics of the phase sum indeed depends on the size of the sampling vol- 
ume explicitly. This feature is naturally explained by the corresponding density structure 
in real space; all the phases are synchronized with each other when a distinctive density 
peak dominates in the given sampling volume, while the phase correlations become weaker 
if several peaks with comparable heights exist. We have confirmed this expectation directly 
by comparing the strength of phase correlations with the ratio of the highest peak to the 
secondary peaks in the sampling volume. Therefore the phase correlation becomes diluted 
as the sampling volume increases because statistically it accommodates a number of peaks 
with comparable heights. 

The fact that the perturbation formula well reproduces the numerical results indicates 
that the information from the phase sum distribution mainly comes from the normalized 
bispectrum amplitude p^ 3 \ Since one can evaluate the bispectrum directly and independently 
from the phase sum distribution, one might wonder if the phase sum statistics is simply 
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Fig. 7. — The relation between and the ratio of the second highest peak density to the 
highest peak density for the 64 cubic data with a box-length of 25/i -1 Mpc (filled circles) 
and the 8 cubic data with a box-length of 50/i _1 Mpc (crosses) based on LCDM cosmological 
model. The ranges of scale \k\\ and |&2| are [0.1,0.2] (Left) and [0.7,0.8] (Right) in units of 
^Nyquist = 27r/(0.78/i- 1 Mpc) and 6 12 = [120°, 130°]. 



Fig. 8. — Same as Figure 6 but for a box with a length of 50h 1 Mpc constituted from the 
eight sub-boxes shown in Figure 6. 
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equivalent to the more conventional bispectrum analysis. Indeed the dependence on the 
sampling volume size may be a key in this context. This implies that any signal of non- 
Gaussianity in the phase correlation disappears in the large-volume limit. We can repeat 
the analysis by systematically decreasing the sampling volume size. Instead, one has to take 
an ensemble average of the phase correlations over many subsamples of a fixed V^mp to 
obtain proper statistical quantities. Otherwise, the phase correlations in a single volume is 
much like the cosmic variance, as the sharpness of the highest peak is weakened in the limit 
V samp — > oo. The situation is similar to the count-in-cells statistics, as the exact number of 
the galaxies in a single cell is attributed to the cosmic variance. Taking an average of the 
weighted number in many cells gives a proper statistics, such as variance, skewness, kurtosis, 
and so forth. The average value and the sample variance with respect to the average as a 
function of V sajnp may give additional information even beyond the bispectrum. 

The distribution of the phase sum may be applied to the real data of the existing galaxy 
catalogues including the 2dF and SDSS, which provides complementary information beyond 
the two-point statistics. Since the phase sum statistics is a new concept in cosmology, we 
still need to explore many other aspects of the statistic that we are likely to have overlooked 
at this point. We are currently working on evaluating the phase-sum distributions of Sloan 
Digital Sky Survey galaxy samples considering all possible systematics by creating simulated 
Mock samples as we did in the previous topological analysis (Hikage et al. 2002, 2003). We 
hope to report the results elsewhere in near future. 

We deeply appreciate I. Kayo and Y. P. Jing for kindly providing a large set of iV-body 
simulation data. We thank K. Yoshikawa for helping us improve our numerical analysis 
routines. C.H. also thanks G. Watanabe, and A. Taruya for useful discussion. Numerical 
computations were carried out at ADAC (the Astronomical Data Analysis Center) of the Na- 
tional Astronomical Observatory, Japan (project ID: yys08a). This research was supported 
in part by the Grant-in-Aid from Monbu-Kagakusho and Japan Society of Promotion of 
Science (12640231, 14102004, 1470157 and 15740151). 
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A. Box-size dependence of the cumulants p( N ) 

According to the analytic formula of the joint PDF of arbitrary sets of the Fourier phase 
9k, the distribution of the phase sum is determined by the cumulants p^ = (0:10:2 • • -o^n)c 
where Oj = fki/y/ (l/fcj 2 ) (A^ > 3 and (• • •) denotes the spatial average). In this section, 
we present the relation of the cumulants p^ to the power spectra and the polyspectra in 
physical unit (independent on the box-size) and explicitly show the dependence of p^ on 
the box-size. First we give the definition of the power spectrum and the polyspectra and 
then apply it to the grid data. 

In a three dimensional density fluctuation field 6(x) = (p(x) — p)/p, the Fourier- 
transform 5(k) of 5(x) in an infinite volume is given by 

f(k) = J d 3 x5(x)e- ik - x . (Al) 

The polyspectra P^ N \ki, ■ ■ •, k^-i) are given by the cumulants 

7(*!i) • • • f(k N )) = (2n) 3 5 3 D (k 1 + • • •fe N )pW(fei, • • ., fe^i), (A2) 



where S D (k) is Dirac's delta function. In this definition, the variance a 2 = {S 2 (x)) is given 
by 

a 2 =[^P(k), (A3) 



2tt 2 

where P(\k\) = p( 2 \k) is the conventional power spectrum. 

In a simulated box with a volume V and a total grid number N, the lattice spacing of x 
is (V/N) 1 / 3 . In a lattice field 5(x) with a discrete set of x's, the discrete Fourier transform 
ffc is computed by 



/ fc = ^E^) e " lfex ' ( A4 ) 

X 

where the lattice spacing in /c-space is 2-k/V 1 ^ 3 . The correspondence between the discrete 
Fourier modes and the continuum Fourier modes is 

lim Wfk = f(k). (A5) 

The cumulants {fk x ■ • • fk N ) c are related to P^ N \ki, ■ ■ •, k N _i) by the relation 

(/ fel • • • h N ) c = V 1 - N ^ l+ ... +kN P^ N \k 1 , ■ ■ ks-i), (A6) 
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where 8% is Kronecker's delta which is 1 (k — 0) and (k ^ 0). The correspondence between 
the Kronecker's delta and the Dirac's delta function is 



lim Vtf = (27r)^(fe). (A7) 

V-*oo,V/N->0 



The cumulants are, therefore, calculated by 

(N) _ ifki ' ' ' /fcjv) 



V(l/fc 1 | z )---(l/ fcjv |^ 

(A8) 



V<I/*I 2 >---<I/*J 2 > 

y/VH-tPikJ ■ ■ ■ P(k N )' 



The box-size dependence of V ' means that higher order terms become dominant 

as the box-size decreases even when k 3 P(k) is small. 
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